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ABSTRACT 

We report the detection of three new exoplanets from Keck Observatory. HD 163607 is a 
metal-rich G5IV star with two planets. The inner planet has an observed orbital period of 
75.29 ± 0.02 days, a semi-amplitude of 51.1 ± 1.4 m s _1 , an eccentricity of 0.73 ± 0.02 and 
a derived minimum mass of Mp sin i= 0.77 ± 0.02 Mj up . This is the largest eccentricity of 
any known planet in a multi-planet system. The argument of periastron passage is 78.7 ± 
2.0°; consequently, the planet's closest approach to its parent star is very near the line of 
sight, leading to a relatively high transit probability of 8%. The outer planet has an orbital 
period of 3.60 ± 0.02 years, an orbital eccentricity of 0.12 ± 0.06 and a semi-amplitude 
of 40.4 ± 1.3 m s _1 . The minimum mass is M P sini= 2.29 ± 0.16 M Jup . HD 164509 is a 
metal-rich G5V star with a planet in an orbital period of 282.4 ± 3.8 days and an eccentricity 
of 0.26 ± 0.14. The semi-amplitude of 14.2 ± 2.7 m s _1 implies a minimum mass of 0.48 
± 0.09 Mj up . The radial velocities of HD 164509 also exhibit a residual linear trend of -5.1 
± 0.7 m s _1 per year, indicating the presence of an additional longer period companion in 
the system. Photometric observations demonstrate that HD 163607 and HD 164509 are 
constant in brightness to sub-millimag levels on their radial velocity periods. This provides 
strong support for planetary reflex motion as the cause of the radial velocity variations. 
Subject headings: planetary systems - stars: individual (HD 163607, HD 164509) 



1. INTRODUCTION 

Over the past 15 years more than 500 extrasolar 
planets have been detected in radial velocity, tran- 
sit, microlensing and most recently, imaging sur- 
veys (http://exoplanet.eu, http://exoplanets.org). 
In step with these discoveries, planet formation 
theories have made use of the observational con- 
straints imposed by the ensemble of exoplanets 
(jFord fc Rasiol[2008T ). Correlations have been un- 
covered between parameters such as the chem- 
ical composition and m ass of host stars and 
the probability of planets (jFischer fc ValentH 120051 : 
I Johnson et all[2007l [20ToT> . 
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Based on the solar nebula model, there was 
an early expectation that planets would reside in 
nearly circular orbits, since eccentric orbits would 
have been quickly dampened through th eir inter- 
actio n with dust in protoplanetary disks (jLissauerl 
1993). Therefore, the detection of significant ec- 
centricity in exoplanet orbits and the unexpected 
parallel between the eccentricity distribution of 
planets and stars continues to be of strong in- 
terest (IMarcyfc Butlerl [19961: ICochran et all [19971: 
IFord fc Rasioll2008[ ) 

One possibility is that non-coplanar orbits 
can be produced by pla net- planet scatter i ng in 
highly eccentric orbit s (fChattcri ee et al.l 120081 : 
Uuric &: Tremainel 120081 ). Among transiting exo- 
planets where the Rossiter-McLaughlin (R-M) effect 
has been measured , approximately on e third seem 
to be misaligned (ITriaud et al.l 120101 : I Winn et"aL1 
l201Ct ISimpson et al.l 120101) . Most of these tran- 
siting systems have low eccentricity orbits; only 
four known transiting systems have eccentricities 
greater than 0.4: H D 17156 (Fischer et al.l[2007 ; 
I Barbieri etall [2007D . HD 80606TNaef et all 12001; 
iFossev et al.l|2009l:IMoutou et al.H2009D. HAT-P-2 b. 
(Bak os et aLll2007l ). and CoR oT-lOb (IBonomo etlbl 
120101) . Using Kepler data, ILissauer et al.l ()2011l ) 
find that nearly coplanar multi-planet systems of 
low mass planets in short-period orbits are quite 
common. However, based on recent observations of 
the fe w misaligned hot- Ju piter systems recently de- 
tected, I Winn et"aLl (|2010D note that misaligned sys- 
tems are preferentially detected around hot stars, 
but the underlying mechanism for misalignment is 
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not clear. 

In this paper we present the detection of two new 
planetary systems orbiting metal rich stars. The 
first is a double-planet system orbiting HD 163607. 
The inner planet, with an orbital period of 75 days 
and an eccentricity of e ~ 0.73, is the most eccen- 
tric planet detected in a multiplanet system to date. 
This eccentricity, combined with the argument of 
periastron passage of 78° , makes the probability of 
transit 8%, 3.7 times higher than it would be if it 
were in a circular orbit with the same orbital period. 
The orbital parameters are precise enough to pro- 
vide good transit ephemeris predictions, and if this 
system transits, then a measurement of the R-M ef- 
fect could shed light on orbital evolution for this 
system. The second system is a 0.5 Mj up planet 
orbiting HD 164509, with an orbital period of 245 
days. 

In section[5]we describe the observations and anal- 
ysis of the two exoplanet systems, followed by a de- 
scription of the photometric followup in section [3l 
In section [4j we describe the properties, observa- 
tions, orbital analysis, transit parameters and pho- 
tometry concerning the host star and two substellar 
companions detected orbiting HD 163607. In sec- 
tion [3] we describe the detection of a jovian mass 
planet orbiting HD 164509. In section [5] we sum- 
marize this paper and discuss the impact of these 
results. 

2. SPECTROSCOPIC OBSERVATIONS AND 
REDUCTIONS 

Doppler measurements for both HD 163607 and 
HD 164509 were obtained at the Keck Observatory 
using the HIRES spectrograph. The intial obser- 
vations were made as part of the N2K program 
(jFischer et al.|[2005T) . which targeted high metallic- 
ity stars to detect short period jovian mass planets. 
Continued observations of these stars are now re- 
vealing longer period and multi-planet systems. 

Radial velocity measurements for the stars on our 
program s are determined with a forward-modeling 
process (jMarcv &: Butler 1992)). The ingredients in 
our model include an intrinsic stellar spectrum (ISS) 
of the star, a high resolution spectrum of the iodine 
cell and a model for the HIRES instrumental pro- 
file (IP). The ISS and the iodine spectrum are mul- 
tiplied together, with the shift of the iodine spec- 
trum as a free parameter to derive the radial veloc- 
ity of the star. This product is then convolved with 
the IP model and binned on the HIRES pixel scale. 
Our Doppler model is broken up into approximately 
seven hundred 2 A chunks from 500 - 600 nm to ac- 
commodate variations in the IP across the detector. 
The single-measurement uncertainty is the weighted 
uncertainty in the mean radial velocities from these 
chunks. 

We derive stellar parameters (Toff, [Fe/H], 
log g and vsmi) using the LTE spectral synthesis 
analysis software Spectroscopy Made Easy (SME ) 
(jValenti fc PiskunovllillM IValenti fc Fischeife OQS) . 
After generating an initial synthe tic model, we it- 
erate between the Y 2 isochrones dDemarque et al. 
2004) and SME model as described by IValenti et al" 
(2009) until agreement in the surface gravity con- 
verges to 0.001 dex. The stellar mass, luminos- 
ity and ages are derived from the Y 2 isochrones 



Table 1 

Stellar Parameters 



HD 163607 HD 164509 



Spectral type 


G5 IV 


G5 V 


V 


8.15 


8.24 


M v 


3.96 


4.64 


B-V 


0.78 


0.66 


BC 


-0.127 


-0.06 


Distance (pc) 


69 (3) 


52 (3) 


T eff (K) 


5543 (44) 


5922 (44) 


log 9 


4.04 (6) 


4.44 (6) 


[Fe/H] 


0.21 (3) 


0.21 (3) 


v sin i(km s *) 


1.49 (50) 


2.4 (5) 


M* (M ) 
R* (Rq) 


1.09 (2) 


1.13 (2) 


1.63 (7) 


1.06 (3) 


L* (L @ ) 


2.3 (2) 


1.15 (13) 


Age (Gyr) 
log R'_ffK 


8.6 (6) 


1.1 (1.0) 


-5.01 


-4.88 




0.164 


0.18 



(jDemarque et al. 2004) and bolometric luminos ity 
corrections are from IVandenBerg &: Cleml (|2003j) . 

Velocity jitter is a combination of astrophysical 
noise and systematic errors that can be misinter- 
preted as dynamical velocities. Starspots are one 
source of stellar jitter. As stars rotate, their spots 
will rise and set over the approaching and receding 
stellar limbs shifting the centroids of spectral lines 
and leading to confusion with d y namic al Doppler 
line shifts. Ilsaacson fc Fischer] (|2010T ) measured 
emission in the Ca II line cores and S*hk values to 
derive log R' HK , the ratio of emission in the core of 
the Ca II lines to the photospheric values. They 
have derived astrophysical jitter measurements as a 
function of B — V color, luminosity class, and excess 
-Shk values and we adopt those stellar jitter mea- 
surements for the stars in this paper. 

To carry out Keplerian modeling of RV and as- 
trometric data, we developed a software package by 
the name of Keplerian Fitting Made Easy (KFME). 
KFME was programmed in IDL and the graphical 
user interface was inspired by th e Systemic Con- 
sole (jMeschiari fc Laughlinl [2010h . The GUI was 
initially developed to fit synthetic radial velocity 
and astrometry data for the Spa ce Interferometry 
Mission (SIM) (ITraub et all 12010) . While SIM has 
been officially discontinued, there is still the pos- 
sibility of using the full potential of KFME, given 
the recent discovery of an exoplanet via astrometry 
(|Muterspaugh et all [20ll ). However, KFME has 
been used here to fit radial velocity data alone. 

KFME displays the data set and allows the user 
to adjust initial parameters for up to seven planets. 
Functions are included for periodogram analysis, 
Levenberg-Marquardt fitting, analysis of false alarm 
probability, and Bootstrap Monte Carlo analysis to 
determine parameter uncertainty. KFME also of- 
fers an automated option that cycles through sev- 
eral values of ui and Tp, retaining the lowest % 2 as 
the best fit solution. 

The false alarm probability (FAP) in KFME is 
calculated with a Monte Carlo simulation. Before 
the MC synthetic data sets are created to calculate 
the FAP for any single planet, any linear trend that 
is present and/or the best-fit Keplerian model for 
additional planets are subtracted from the veloci- 
ties. Then, for each MC trial, the observation times 
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are kept fixed, while the associated velocities are 
scrambled (replacement, with redraw of the same 
value is allowed) . Each MC trial of scrambled veloc- 
ities is then blindly fit with a Levenberg-Marquardt 
least-squares minimization algorithm with a Kep- 
lerian model. The lowest xt f° r each test case is 
stored in an array. After N trials (where N is set 
by the desired precision in the FAP) the array of 
xt values is sorted and a comparison with the orig- 
inal xt °f the unscrambled velocities is made to 
this ranked set. The FAP assesses how frequently 
the xt °f the scrambled velocities is lower than the 
xt °f the original, unscrambled velocities. For ex- 
ample, if the scrambled data yielded periodograms 
with peaks of comparable height as the original peak 
in 10 out of 100 trials, the FAP of 10/100 = 0.1 
reflects the fact that spurious signals could have oc- 
curred 10% of the time. 

For multiplanet systems, signals can be individu- 
ally displayed. KFME allows fitting for systematic 
velocity offsets (e.g., different observatories, differ- 
ent detectors, etc.) and the inclusion of jitter, added 
in quadrature to the formal errors. After fitting for 
one or more planets, residual velocities can be dis- 
played and analyzed. 

One advangage of KFME is that it engages 
the human brain, which is often good at dis- 
cerning global patterns, to approximate the ini- 
tial conditions for generating a Keplerian model 
with a Levenberg-Marquardt fitting algorithm. The 
Levenberg-Marquardt algorithm then polishes this 
approximate solution into a low xt fit. The 
best fit can then be used to generate transit pa- 
rameters (ingress, egress, duration, time of cen- 
ter transit, and probability). KFME also col- 
lects all of our orbital modeling tools into one 
package. KFME can be downloaded at exoplan- 
ets.astro.yale.edu/KFME, and requires IDL in or- 
der to run. 

3. PHOTOMETRIC OBSERVATIONS AND 
REDUCTIONS 

We acquired photometric observations of 
HD 163607 and HD 164509 with the T12 0.80 m 
automatic photometric telescope (APT) at Fairborn 
Observatory. The T12 APT and its two-channel 
photometer measure photon count rates simulta- 
neously through Stromgren b and y filters. T12 
is essentiall y identical t o the T8 0.80 m APT 
described in lHenryl (|1999l ). 

The two program (P g ) stars HD 163607 and 
HD 164509 were each observed differentially with 
respect to two nearby comparison stars (CI and 
C2). The two comparison stars for HD 163607 
were HD 169352 (V = 8.01, B -V = 0.44, F2) 
and HD 165700 (V = 7.79, B - V = 0.46, F8); 
comparison stars for HD 164509 were HD 166073 
(V = 6.99, B — V = 0.45, F7 IV) and HD 165146 
(V = 7.57, B — V = 0.45, F0). The differential 
magnitudes P g — CI, P g — C2, and C2 — CI were 
computed from each set of differential measures. 
The observations were corrected for extinction and 
transformed to the Stromgren photometric system. 
To improve the precision of our brightness measure- 
ments, we averaged the b and y differential magni- 
tudes into a single (b + y)/2 "passband", which we 



Table 2 

Radial Velocities for HD 163607 



JD 


RV 




Shk 


-2440000 


(m s — *) 


(m s~ 




13570.8646 


-23.74 


1.18 


0.164 


13575.9339 


-13.41 


1.24 


0.161 


13576.8468 


-7.42 


1.30 





14247.9039 


38.39 


1.13 


0.164 


14249.8472 


46.49 


1.27 


0.164 


14250.9825 


47.65 


1.18 


0.165 


14251.9291 


52.15 


1.13 


0.164 


14285.8917 


-7.61 


0.97 


0.164 


14318.9150 


23.03 


0.83 


0.164 


14339.8773 


-44.90 


1.05 


0.164 


14343.8563 


-41.10 


0.93 


0.165 


14345.7946 


-31.49 


1.37 


0.167 


14345.8208 


-33.73 


1.05 


0.164 


14549.0883 


-6.88 


1.18 


0.164 


14634.9349 


14.16 


1.00 


0.164 


14639.0081 


-84.13 


1.00 


0.168 


14641.9689 


-90.54 


1.12 


0.164 


14674.8896 


-49.94 


1.02 


0.165 


14689.9213 


-36.62 


1.12 


0.165 


15016.0415 


-79.07 


0.96 


0.164 


15041.9308 


-39.20 


1.14 


0.165 


15044.0017 


-42.11 


1.08 


0.164 


15073.7742 


-0.64 


1.11 


0.165 


15082.7333 


26.03 


1.09 


0.163 


15111.7501 


-35.94 


1.20 


0.165 


15135.7005 


-2.87 


1.13 


0.164 


15163.6885 


17.28 


1.30 


0.161 


15172.6923 


-42.17 


1.19 


0.162 


15229.1489 


46.57 


1.06 


0.163 


15232.1671 


56.60 


1.14 


0.161 


15256.1669 


-7.61 


1.16 


0.162 


15261.0850 


-1.71 


1.13 


0.164 


15285.0537 


28.62 


1.13 


0.154 


15313.9685 


50.43 


1.13 


0.163 


15319.0701 


-17.97 


1.22 


0.161 


15351.1269 


11.69 


1.02 


0.163 


15373.8050 


43.04 


1.09 


0.164 


15399.9669 


-12.62 


1.02 


0.163 


15436.7391 


27.49 


1.07 


0.157 


15455.7818 


53.65 


1.00 


0.164 


15486.7739 


1.75 


1.04 


0.164 


15500.6979 


10.70 


1.16 


0.164 


15521.7059 


32.83 


1.04 


0.163 


15606.1749 


40.56 


1.09 


0.164 


15606.1796 


41.72 


1.11 


0.163 


15613.1249 


61.89 


1.84 


0.158 


15634.0875 


-17.56 


1.03 


0.164 


15637.0856 


-17.19 


1.06 


0.163 


15669.0353 


14.99 


1.07 


0.164 


15671.0957 


18.70 


1.08 


0.164 


15700.8908 


-44.97 


1.15 


0.164 


15728.9113 


-16.27 


1.16 


0.164 


15763.7659 


40.29 


1.00 


0.164 



designate by. Typical precision of a single by obser- 
vation is ~ 0.0015 — 0.0020 mag, a s measured for 
pairs of constant stars. IHenrvl (| 1999T ) provides addi- 
tional details on the operation of the APT, observ- 
ing and data reduction procedures, and precision of 

t he data. 

IQueloz et all (f200l and IPaulson et all (f2004l ) 
have demonstrated how rotational modulation of 
starspots on active stars can result in periodic ra- 
dial velocity variations that mimic the presence of 
a planetary companion. Thus, the precise APT 
brightness measurements are valuable for distin- 
guishing between activity-related RV changes and 
true reflex motion of a star caused by a planetary 
companion. 

4. HD 163607 
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HD 163607 (HIP 87601) is a G5 subgiant at a 
distance of 69 ± 3 pc c alculated fr om the Hipparcos 
par allax measurement j ESA 1997]) and revised cata- 
log ([van Le cuwen 2008) . We adopted the Hipparcos 
V-band magnitude and color of V = 8.15 and B — V 
= 0.78. With a bolometric correction of -0.127, this 
gives the absolute visual magnitude of My = 3.96. 

An iodine-free "template" spectrum of HD 163607 
was analyzed by iterating SME models with Y 2 
isochrones to derive the following stellar parame- 
ters: T eS = 5543 ± 44 K, [Fe/H] = 0.21 ± 0.03 dex, 
projected stellar rotational velocity, vsini— 1.5 ± 
0.5 km s , \ogg= 4.04 ± 0.06. The isochrone anal- 
ysis also yielded an age of 8.6 ± 0.6 Gyr, a stellar 
radius 1.63 ± 0.07 R© and a luminosity of 2.3 ± 
0.2 L Q . HD 163607 has low chromospheric activ- 
ity with logR' HK = —5.01 and an estimated stellar 
jitter of 2.6 m s _1 . We do not see a correlation be- 
tween activity and the measured radial velocities. 
The stellar properties for HD 163607 are summa- 
rized in Table [TJ 



The uncertainties quoted in this section and listed 
in Table [JJ come from 10 3 realizations of this Boot- 
strap Monte Carlo routine. 
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Figure 1. Phase-folded radial velocities for HD 163607b, 
the inner planetary companion in this system, which has an 
orbital period of 75.29 days. The theoretical velocities from 
the outer planet have been subtracted from the observations 
and the curve shown is the best-fit model for this system. 



4.1. Orbital Solution 

Observations of HD 163607 began in July of 2005. 
The 51 observations of this star are listed in Table [1] 
and have a median velocity precision of 1.11 m s . 
We analyzed the radial velocity data for HD 163607 
using KFME and included jitter of 2.6 m s _1 . The 
single planet model had a period of 1500 days, how- 
ever, the residuals to this fit had an rms of 31 m s _1 
and significant power at 75 days. The addition of 
a second planet in the Keplerian model reduced the 
xt fit to 1.03 with an rms of 2.9 m s _1 . Once the 
best-fit two planet Keplerian model was attained us- 
ing KFME, a Bootstrap Monte Carlo routine that 
is also built into KFME was used to derive un- 
certainties for the orbital parameters. The Boot- 
strap Monte Carlo routine, similar to the routine 
employed to calculate the FAP described in section 
[21 subtracts the best-fit model from the velocities, 
scrambles the residuals and their associated uncer- 
tainties, adds them back to the model and refits. 



Figure 2. Radial velocities for HD 163607 with the theo- 
retical velocities for the inner planet removed from the ob- 
servations. The curve shown is the best-fit Keplerian model 
for HD 163607c. 

The inner planet has a best-fit period, P, of 
75.29 ± 0.02 days and velocity amplitude, K, of 
51.1 ± 1.4 m s _1 . The best-fit eccentricity, e, for 
this planet is 0.73 ± 0.02. Assuming a stellar mass 
of 1.09 Mq, the derived minimum planet mass is 
M P sin i= 0.77 ± 0.04 M Jup . The high eccentric- 
ity and favorable argument of periastron passage cj, 
conspire to increase the probability of transit to 8%; 
much larger than it would have been if the inner 
planet was in a circular orbit. The observed radial 
velocities for the inner planet as a function of or- 
bital phase are shown in Figure [TJ In this figure the 
Keplerian model for the outer planet has been re- 
moved and the Keplerian model for the inner planet 
is plotted as the solid line. 

The outer planet, HD 163607c, has a best-fit pe- 
riod, P, of 1314 ± 8 days, a velocity semi-amplitude 
of 40.4 ± 1.3 m s _1 and an eccentricity, e, of 0.12 
± 0.06. The derived mass is 2.29 ± 0.16 M Jup and 
the semi-major axis of the orbit is 2.42 ± 0.01 AU. 
In Figure [2 the Keplerian model fit to the data 
for HD 163607c is shown with the inner planet sub- 
tracted from the Doppler measurements. MCMC 
simulations to confirm our calculation of the best-fit 
orbita l paramete r s were carried out using the meth- 
ods of Ho u et aTl (|2011[) . and the resulting posterior 
PDFs for several orbital parameters can be seen in 
Figure [3J There was excellent agreement between 
the MCMC results and the frequentist approach us- 
ing KFME. The orbital parameters of both planets 
obtained using KFME are summarized in Table [TJ 

Similar to the scrambled velocities method de- 
scribed in Marcv et alj (|2005l ). the false alarm prob- 
abilities (FAPs) for each planet were estimated by 
creating 10 4 synthetic data sets by drawing, with 
replacement, velocities and their associated errors 
from the data, and placing these at the actual obser- 
vation times. A Levenberg-Marquardt least-squares 
minimization to a Keplerian model was then per- 
formed for each synthetic data set, and the xt dis- 
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Table 3 

Orbital parameters for the three planets described in this 
work. 



Parameter 


HD 163607b 


HD 163607c 


HD 164509b 


p(d) 


75.29 ± 0.02 


1314 ± 8 


282.4 ± 3.8 


T P (HJD*) 


14185.00 ± 0.24 


15085 ± 880 


15703 ± 30 


t c (HJD*) 


15841.59 ± 0.24 


17074 ± 15 


15498 ± 22 


e 


0.73 ± 0.02 


0.12 ± 0.06 


0.26 ± 0.14 


U) 


78.7 ± 2.0 


265 ± 93 


324 ±110 


K (m s" 1 ) 


51.1 ± 1.4 


40.4 ± 1.3 


14.2 ± 2.7 


a (AU) 


0.36 ± 0.01 


2.42 ± 0.01 


0.875 ± 0.008 


M sini (Mj U p) 


0.77 ± 0.04 


2.29 ± 0.16 


0.48 ± 0.09 


7 (m s _i ) 


-15.7 ± 0.5 




8.9 ± 2.1 


dvdt (m s yr ) 







-5.1 ± 0.7 


N obs 


51 




41 


Jitter (m s ) 


2.6 




3.2 


rms (m s _1 ) 


2.9 




4.9 


X v 


1.03 




2.04 



HJD - 2,440,000 




16.1 75.1S 75.1* 75.25 75. S 75 .35 75 .4 



4 ^^K^ EGcentrldty, 

0?58 0.7 0.72 7a 07 5 0.78 0.8 




0-OS 0.1 0.15 0.2 0.25 



Figure 3. MCMC histograms showing the posterior PDFs 
for several orbital parameters for HD 163607b (left) and 
HD 163607c (right). The very narrow PDFs for HD 163607b 
reflect the excellent phase coverage near periastron passage 
as seen in Figure [l] 

tributions for both planets orbiting HD 163607 are 
shown in Figure 01 The \t f° r the unscrambled 
data set is indicated by the downward pointing ar- 
row in each plot, which is much lower than the xt °f 
any of the 10 4 scrambled velocities for both planets 
in this system. This results in a FAP of < 0.01%, in- 
dicating that the observed signal is not due to noise, 
but rather a coherent periodicity that is well-fit with 
a Keplerian model. 



4.2. Transit Prediction 

The proximity of hot Jupiters to their host stars 
results in an increased transit probability compared 
to longer period planets. Although the orbital pe- 
riod of HD 163607b is 75 days, the high orbital ec- 
centricity moves the planet close to the star during 
periastron passage. Since the periastron passage of 
HD 163607b is serendipitously oriented very close 
to our line of site, the probability of transit is much 



larger than it would be if the planet was in a circular 
orbit. To illustrate this point, the orbital configu- 
ration for HD 163607b is shown in Figure [5] The 
calculation of the time of center transit, ingress and 
egress are particularly helpful for planning photo- 
metric follow-up for longer period planets. Below 
we describe a technique for calculating parameters 
to aid in photometric followup. 

4.2.1. Transit Center, Probability and Duration 

The term time of center transit, t c , as defined here 
is the time when the planet is nearest the center 
of the disk of its parent star as seen from Earth. 
We calculated t c by using the best-fit orbital pa- 
rameters to produce an array of Cartesian coor- 
dinates, giving the position of the planet relative 
to the host star in the inclination-projected orbital 
plane. The inclination-projected orbital plane was 
then rotated into the plane defined by the Earth, 
ascending, and descending nodes using the Thiele- 
Innes coordinates. The ingress and egress were then 
derived using the coordinates where the planet en- 
tered and exited the disk of the parent star. The 
true anomaly corresponding to each position can 
then be found and in turn the mean anomaly calcu- 
lated. The mean anomaly, M, can then be used to 
find the time since periastron passage using Equa- 
tion □ 



PM 
~2^~ 



(1) 



where P is the orbital period, M is the mean 
anomaly and t is the time since periastron passage. 
Finding the difference in the mean anomalies for 
egress and ingress gives the transit duration, td- The 
transit center is then given by equation [2] 



tr 



t + Tp 



(2) 



Analytically, the time of central transit fol- 
lows from the definition of the true anomaly 
and the argument of pe r iastron, as d escribed by 
Tingl ev fc Sacketti (f2005h : iKanel l|200l : IKane et al.1 
(2009). The probability that the planetary com- 
panion will transit is given by ICharbonneau et "aLl 
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Figure 4. To determine the false alarm probability (FAP) for each planet, the velocities were scrambled 10,000 times and 
each set of velocities was blindly fit with a Keplerian model. The above figures show the distributions of these realizations as 
a function of \l for HD 163607b (left), HD 163607c (middle) and HD 164509b (right). The dashed arrow shows the xl for the 
unscrambled velocities, which is well to the left of each distribution. The FAP for each planet was < 0.01%. 



Table 4 

SUMMARY OF PHOTOMETRIC OBSERVATIONS FROM THE T12 APT 



Program 
Star 
(1) 


Date Range 
(HJD - 2.400,000) 
(2) 


Duration 
(days) 
(3) 


N ob3 
(4) 


a(C2 - Cl) by 
(mag) 
(5) 


a(P g - ClC2) by 
(mag) 
(6) 


Planet 
(7) 


Orbital Period 

(days) 
(8) 


Scmiamplitude 
(mag) 
(9) 


HD 163607 


54381-55499 


1118 


156 


0.0019 


0.0015 


b 


75.29 


0.0002 ± 0.0002 














c 


1314 


0.0006 ± 0.0003 


HD 164509 


54377-55479 


1102 


168 


0.0018 


0.0018 


b 


282.4 


0.0006 ± 0.0002 



3 
< 




0.2 



0.0 



-0.2 -0.4 
Position (AU) 



-0.6 



Figure 5. A view of the HD 163607b planetary orbit. The 
direction of the observer is indicated in the bottom of the 
illustration. The dots indicate the ingress and egress for the 
primary and secondary transits. 



(|2007f ) and is reproduced in Equation [3] The last 
factor in square brackets in this equation shows that 
HD 163607 is 3.7 times more likely to transit its par- 
ent star than if it were in a circular orbit with the 
same semi-major axis. 



- Circular Orbit 
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o Transit 
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1 10 100 1000 

Orbital period (days) 



10000 



Figure 6. The probability of transit versus orbital period 
for the known distribution of exoplanets. The solid line is 
the probability for a circular orbit of a Jupiter radius planet 
around a one solar radius star. The filled circles represent 
planets detected through the radial velocity method that 
have not been observed to transit. The open circles rep- 
resent planets that are known to transit and the star symbol 
represents HD 163607b. 



p = 4.5*10 



-3 f lAU 
V a 



+ R 



pi 



Rc. 



1-e 2 



(?) 

In Equation[3J a is the semi-major axis, i?* is the 
radius of the star, R p i is the radius of the planet, e is 
the eccentricity of the orbit and u> is the argument 
of periastron passage, which is the angle between 
the ascending node and the position of the planet at 
periastron. Figure [6] shows the probability of transit 
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for the known distribution of planets when taking 
the eccentricity and argument of periastron passage 
into account. The solid line shows the probability of 
transit as a function of orbital period for a Jupiter 
radius exoplanet transiting a solar radius star in a 
circular orbit. The filled circles are planets that 
have only been detected through the radial velocity 
method and the open circles are planets that have 
been observed to transit their parent stars. For both 
the filled and unfilled circles, the parent star's radius 
is set to solar and each planet's radius is set to the 
radius of Jupiter. The blue star shows where the 
inner companion orbiting HD 163607 falls on this 
plot, illustrating how the favorable orientation and 
high eccentricity enhance the probability of transit. 

Using the equations described in this section we 
calculated the next time of center transit to be t c 
= 2455841.59 ±0.24 UT, the duration is roughly 
5 hours, the depth is ~ 0.4% (~4 mmag) and the 
probability of transit is 8% for HD 163607b. The 
uncertainty in t c is determined while calculating 
the uncertainty in the orbital parameters. A boot- 
strap monte carlo routine is employed that subtracts 
the best fit Keplerian model, scrambles the residu- 
als and their associated errors, adds these residuals 
back to the theoretical curve and refits. The new 
best fit orbital parameters are then used to deter- 
mine t c . The standard deviation of the resulting t c 
array after 10 3 realizations is what is stated for the 
uncertainty in t c . The predicted transit window is 
defined in this work as the time between ingress and 
egress plus or minus one sigma. None of our radial 
velocity spectroscopic observations happened to fall 
within the predicted transit window. 

4.3. Photometry 

Photometric results for HD 163607 are given in 
the first two rows of Table [TJ The observations 
were acquired between 2007 October 3 and 2010 
October 29. Column 5 gives the standard devia- 
tion of the comparison star differential magnitudes 
C2 — CI as 0.0019 mag. This is typical for constant 
stars with this telescope. Periodogram analysis of 
the comparison star observations did not detect any 
significant periodicity between one and 200 days, so 
both comparison stars are constant to the level of 
the photometric precision. To improve the preci- 
sion of the HD 163607 observations, we computed 
the P g — CI and P g — C2 differential magnitudes 
in the (b + y)/2 passbands and then took the mean 
of those two differential magnitudes. This results 
in differential magnitudes in the sense HD 163607 
minus the mean brightness of the two comparison 
stars, which we designate as P g — C1C2. The stan- 
dard deviation of the P g — C1C2 observations is 
given in column 6 of Table [T] The small value of 
0.0015 mag and the absence of any significant peri- 
odicity in the P g — C1C2 observations indicates that 
HD 163607 is also constant to the limit of precision. 

We computed a least-squares sine fit of the 
P g — C1C2 observations phased separately to the 
75.29- and 1314-day periods of planets b and c, re- 
spectively. Planet b exhibits a semiamplitude of 
0.0002±0.0002 mag, which is consistent with zero to 
high precision. Planet c has a slightly larger semi- 
amplitude of 0.0006 ± 0.0003 mag. We note that 



Table 5 

Radial Velocities for HD 164509 



i n 


T-iV 
rtv 




JHK 


o a a r\r\r\r\ 

-Z44UUUU 


(m s ) 


/, „ — 1 \ 

(m s ) 




13570.8176 


0.87 


1.47 


0.187 


13576 0255 


16.38 


1.56 


0.180 


13576.8579 


14.03 


1.43 




14251.0462 


4.61 


1.28 


0.179 


14318 8554 


27.68 


1.14 


0.182 


14339.7544 


20.51 


1.25 


0.184 


14343.8855 


21.89 


1.16 


0.183 


14633.9129 


12.49 


1.33 


0.183 


14634.8953 


12.71 


1.37 


0.184 


14635.9216 


10.06 


1.28 


0.183 


14637 0538 


3.71 


1.36 


0.183 


14637.9282 


11.22 


1.22 


0.182 


14638.9721 


5.64 


1.29 


0.182 


14639.9825 


4.53 


1.18 


0.180 


14641 9386 


4.27 


1.36 


0.179 


14644.0757 


3.78 


1.45 


0.178 


14674.8405 


4.47 


1.29 


0.183 


14688 8406 


-7.39 


1.36 


0.182 


14930.0681 


7.13 


1.38 


0.192 


14930.0700 


5.49 


2.16 


0.188 


14956.1303 


-0.07 


1.27 


0.189 


14964.0480 


6.96 


1.23 


0.189 


14985.1079 


-1.27 


1.47 


0.190 


15016.9568 


-16.55 


1.27 


0.188 


15019.0135 


-18.16 


1.50 


0.186 


15041.9070 


-10.20 


1.33 


0.187 


15043.7946 


-11.64 


1.37 


0.187 


15048.7836 


-14.61 


1.41 


0.188 


15077.7308 


-11.39 


1.18 


0.182 


15106.7446 


2.18 


1.29 


0.183 


15135.7359 


8.32 


1.33 


0.185 


15229.1718 


-1.50 


1.32 


0.186 


15232.1564 


-9.79 


1.29 


0.180 


15256.1429 


-17.10 


1.25 


0.172 


15286.1172 


-14.85 


1.29 


0.185 


15350.9342 


-15.92 


1.19 


0.185 


15380.7878 


-13.97 


1.37 


0.179 


15455.7291 


3.08 


1.21 


0.182 


15486.7349 


-4.27 


1.30 


0.181 


15500.7174 


1.43 


1.28 


0.182 


15521.6869 


-12.83 


1.46 


0.178 


15522.6823 


-15.26 


1.42 


0.178 


15607.1406 


-24.51 


1.25 


0.183 


15636.1381 


-9.79 


2.91 


0.179 


15668.0865 


-17.46 


1.14 


0.178 


15669.0153 


-14.03 


1.27 


0.176 


15671.0647 


-7.97 


1.34 


0.177 


15701.1294 


3.91 


1.37 


0.178 


15723.9194 


3.25 


1.42 





our observations do not quite cover one full cycle 
of the 1314-day RV variations. However, given the 
very low standard deviation of the P g — C1C2 ob- 
servations, we conclude that HD 163607 is also con- 
stant to high precision on the period of planet c. 
These very tight limits of brightness variability on 
the 75.29- and 1314-day RV periods strongly sup- 
port the interpretation that the RV variations in 
HD 163607 are due to stellar reflex motion in grav- 
itational response to planets b and c. Our photo- 
metric observations are too few for transit searches 
for the two planets. 

5. HD 164509 

HD 164509 (HIP 88268) is a G5 main sequence 
star at a distance of 52 ± 3 pc calc ulated from 
the Hipparcos par allax measurement {ESA 1997; 
Ivan Leeuw cn 2008). We adopted the V-band mag- 
nitude and color from the revised Hipparcos cat- 
alog of V = 8.24 and B - V of 0.66, and derive 



Giguere et al. 



I 10 r / 




-30 1= 

2004 2005 



2006 2007 2009 
Time of Observation 



2010 2011 



Figure 7. The precise Dopplcr velocities for HD 164509 are 
shown here as the black circular points with associated 1 a 
uncertainties. The solid curve is the theoretical Keplerian 
model that best fits this system. 

an absolute visual magnitude of My = 4.64, which 
includes a bolometric correction of -0.06. Spectro- 
scopic analysis of HD 164509 yields: T eS = 5922 ± 
44 K, [Fe/H] = 0.21 ± 0.03 dex, vsini= 2.4 ± 0.5 
km s _1 and \ogg= AAA ± 0.06. The isochrone anal- 
ysis yields a mass of 1.13 ± 0.02 M Q , an age of 1.1 
± 1 Gyr, a stellar radius of 1.06 ± 0.03 R Q and a 
luminosity of 1.15 ± 0.13 L Q . The star has modest 
chromospheric activity, with logi?^ = —4.88 and 
llsaacson fc Fisc her (2010) estimate a stellar jitter of 
3.2 m s . The stellar properties described above 
for HD 164509 are summarized in Table [T] 

5.1. Orbital Solution 

Observations of HD 164509 began in July of 2005 
at Keck Observatory with the HIRES spectrometer; 
and 41 observations of this star now span a time 
baseline of 5 years. The observations are listed in 
Tableland the median velocity error for HD 164509 
is 1.32 m s _1 . We carried out Keplerian modeling 
for HD 164509 using KFME after adding jitter of 
3.2 m s _1 in quadrature with the formal velocity 
errors. 

Our best-fit model is for a single planet with a pe- 
riod of 282.4 ± 3.8 days, orbital eccentricity of 0.26 
± 0.14, a radial velocity amplitude of 14.2 ± 2.7 
m s _1 and a linear trend of -5.1 ± 0.7 m s _1 yr -1 . 
With the assumed stellar mass of 1.13 ± 0.02 Mq, 
we derive a mass for the planet of Mp sini= 0.48 ± 
0.09 Mj up . The rms to our model fit is 4.9 m s _1 ; 

with the assumed jitter, the xl f° r this fit is 2.04, 
suggesting that either the jitter has been underesti- 
mated or there are additional weak dynamical sig- 
nals that have not been adequately modeled with 
a single planet fit. Figure [7] shows the best fit Ke- 
plerian model for the data and the orbital param- 
eters are summarized in Table [T] Similar to the 
HD 163607 syst em, we car r ied ou t MCMC simu- 
lations following IHou et all (|2011[ ). The resulting 
posterior PDFs can be seen in Figure [51 which were 
again in excellent agreement with the best fit so- 
lution from KFME. Just as with HD 163607, an 
FAP analysis was carried out for HD 164509b. Out 
of the 10 synthetic data sets created, not a single 



1Sj 



10 



x 10 




10 12 14 16 18 2D 22 24 



^10 




Figure 8. MCMC posterior PDFs for several orbital pa- 
rameters for HD 164509b. 

set had a lower xt fit than the unscrambled veloci- 
ties, leading to an FAP of < 0.01%. The right-most 
plot in Figure [4] shows the xt distribution for the 
10 4 synthetic data sets created to estimate the FAP 
of HD 164509b, and the downward pointing arrow 
shows the xt f° r the best fit to the unscrambled 
velocities. 

5.2. Photometric Observations 

Results for HD 164509 are given in the third row 
of Table [TJ These observations were acquired be- 
tween 2007 October 26 and 2010 October 9. The 
standard deviation of the comparison star observa- 
tions is 0.0018 mag, consistent for constant stars. 
Periodogram analysis of these C2 — CI compari- 
son star observations finds no significant periodic- 
ity. We compute the P g — C1C2 differential mag- 
nitudes as we did for HD 163607 above. The re- 
sulting low standard deviation of 0.0018 mag, along 
with the absence of any significant periodicity in 
the P g — C1C2 observations, demonstrates that 
HD 164509 is constant to high precision. 

A least-squares sine fit of the P„ — C1C2 obser- 
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vations phased to the 282.4-day RV period gives a 
semi-amplitude of 0.0006 ± 0.0002 mag. This tight 
limit to brightness variability further supports the 
existence of HD 164509b. Again, our observations 
are too few for any transit search. 

6. SUMMARY & DISCUSSION 

In this paper, we present Doppler velocities 
for two metal rich stars from the N2K survey: 
HD 163607 and HD 164509. The velocities for the 
G5 IV star, HD 163607 are well-fit with a two- 
planet Keplerian model with xt = 1-03 and resid- 
ual rms of 2.9 m s^ 1 . Based on the 75 and 1314 
day orbital periods and semi-amplitudes listed in 
Table [TJ the derived minimum masses are Mp sini 
= 0.77 M Jup and M P sinz = 2 29 Mj up , respec- 
tively. Ilsaacson fc Fischerl (|2010[ ) estimate the jit- 
ter for a subgiant star of this color to be 4.3 m s — 1 , 
which gives a xl 01 0.4. However, Hipparcos data 
show that this star is only mo destly evolved off of 
the m ain sequence. Using the Isa acson fe Fischerl 
(|2010D jitter for a main sequence star of this color 
of 2.6 m s -1 gives a xt of ~ 1, indicating this lower 
value is a much better estimate of the uncertainty, 
which is why we chose to use that lower estimate 
for the uncertainty for this work. 

The velocities for HD 164509 are best fit by a 
single Keplerian model with a period of 282 days 
and a linear trend. The planet is in an orbit with 
modest eccentricity, e = 0.26, and has a minimum 
mass, Mp sini— 0.48 Mj up . However, the Keplerian 
model has a RMS of 4.9 m s" 1 and xl of 2.04, 
suggesting that the model does not fully describe 
our data. In this case, the stellar jitter may have 
underestimated astrophysical noise, or an additional 
companion is contributing to the residual velocities. 

The HD 163607 system is particularly interesting 
for a number of reasons. The inner planet, with 
an eccentricity of 0.73, is the most eccentric planet 
in a multiplanet system. The average eccentricity of 
plane ts in multiplanet systems is 0.22 ( Wrig ht et al.1 
2009). It is not clear why this system harbors such 
a high eccentricity inner planet. Could the inner 
planet have been scattered inwards by the outer 
planet? Or possibly a third planet that was ejected 
from the system? Is the Kozai mechanism respon- 
sible, due t o the outer planet or a distant stellar 
companion (|Nasasawa et all [20081? Are we seeing 
a high-eccentricity snapshot of a system that os- 
cillates by dynamical interactions between low and 
high eccentricity states? 

Because of the high eccentricity (e = 0.73) and 
orbital configuration with respect to the Earth (w 
= 79°) of HD 163607b, the probability of transit 
for the inner planet is much higher than it would 
have been if the orbit was circular. Planet-planet 
scattering simulations suggest that highly eccentric 
systems might also have large spin-orbit misalign- 
ments, yet there are only a handf ul of known sys- 
tems to test this theoretical re sult ([Chattcricc et al.l 
l200llJuric fc Tremainell2008D . If HD 163607b were 
to transit, spectroscopic observations taken during 
transit for this bright and highly eccentric system 
would allo w for the calculatio n of the spin-orbit 
alignment (IQueloz et a ll [20011: IQhta et all [20051: 
I Winn et "all 120091 ) . making this an excellent sys- 



tem to test theoretical predictions. While we have 
attempted to detect a transit event for the inner 
planet orbiting HD 163607, the combination of the 
orbital period being relatively long and the orbital 
period being a non-integer multiple of a day have 
made these attempts quite difficult. We encour- 
age members of the community to search for transit 
events of this inner planet. 
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